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Abstract In this paper we present a method for deriving effective one-dimensional models based on the 
matrix product state formalism. It exploits translational invariance to work directly in the thermodynamic 
limit. We show, how a representation of the creation operator of single quasi-particles in both real and 
momentum space can be extracted from the dispersion calculation. The method is tested for the analytically 
solvable Ising model in a transverse magnetic field. Properties of the matrix product representation of the 
creation operator are discussed and validated by calculating the one-particle contribution to the spectral 
weight. Results are also given for the ground state energy and the dispersion. 

PACS. 02.70.-c Computational techniques; simulations - 75.10.P Spin chain models - 05.10.Cc Renor¬ 
malization group methods 


1 Introduction 

Strongly correlated quantum systems remain a great chal¬ 
lenge in condensed matter theory. In many cases, we have 
to rely on numerical tools to make predictions from a mi¬ 
croscopic model that can be compared to experiment. Un¬ 
fortunately the exponential growth of the Hilbert space 
dimension severely limits the size of systems that can be 
treated exactly. 

However, to gain understanding of the physics in quan¬ 
tum systems, often the properties of the ground state and 
a few excited states go a long way. Therefore, a variety of 
tools have been developed to separate the low energy part 
of the Hilbert space from the rest. In one dimension, this 
is most notably the density matrix renormalization group 
(DMRG) [1], which is a variational ansatz. There are vari¬ 
ous extensions of the original DMRG, many of which work 
within the framework of matrix product states (MPS) [2] 
or extensions thereof [3,4]. 

Another approach that extends more readily to higher 
dimensions is the semi-numerical method of continuous 
unitary transformations (CUTs) [5-9]. The idea behind 
these CUTs is to partially diagonalize the Hamiltonian 
and derive an effective model (explained below) in sec¬ 
ond quantization for the low-energy sector of the Hilbert 
space. This effective model can then be used to calculate 
physical properties of the system. CUTs have been applied 
successfully in a variety of cases [10-12]. 

In our view, it is a promising long-term goal to es¬ 
tablish effective models for strongly correlated systems 
in terms of the elementary excitations with the ground 
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state being the vacuum. Assuming translational invariance 
for infinitely large systems, i.e., systems in the thermody¬ 
namic limit, the momentum space notation is convenient 
due to momentum conservation and mutual orthogonal¬ 
ity of momentum eigenstates with different momenta. The 
goal of our approach is to map a microscopic lattice model 
such as the transverse held Ising model (TFIM), see Eq. 
(2), to an effective Hamiltonian for the low-energy physics 
in the quasi-particle picture 

n cS = e 0 + ^2 a h a ?i ( la ) 

<?1 

+ ~7f 53 [-^92,93 a 92 a L a 9i + h- c -] <5(32+93,91 (lb) 

V ^ 91,92,93 

+ £ 53 Vqz'qi \ a \ 1 a \ 3 a q2 a qt\ ^qi+q2, 93+94 (1°) 

<71,92, 

93,94 

+ [higher terms] . 

In T~L cS the second quantization operator d' q . creates a 
quasi-particle with momentum qi and a qi annihilates one. 
The ground state energy is labeled E 0 and the disper¬ 
sion relation uj q . The and V// 1 // 2 are the matrix 

elements for quasi-particle decay and two-particle inter¬ 
action, respectively. Momentum conservation is included 
through the Kronecker S symbols. At this point we do not 
include further algebraic properties, since these depend on 
the model under consideration. 

By the effective model (1), the fundamental physical 
properties of complicated microscopic models can be re¬ 
duced to a level which makes further quantitative studies 
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possible. The above mentioned continuous unitary trans¬ 
formations (CUTs) represent a systematic tool which suc¬ 
cessfully achieves this aim [6,13]. They work particularly 
well if the chosen starting point, the so-called reference 
state, is already close to the true ground state of the sys¬ 
tem. But if this is not the case or if it is even not possible 
to find a technically tractable starting point with the rel¬ 
evant symmetries the CUTs cannot be applied efficiently. 
In such situations a completely numerical approach is ap¬ 
pealing because it can tackle a problem at hand in an 
unbiased fashion. 

This is where the MPS representation comes into play. 
It is extremely efficient in finding ground states and cer¬ 
tain excited states. We will show how one can use MPS 
to numerically define and derive local creation and anni¬ 
hilation operators at site i. Thereby, we provide the 
first steps towards representations such as (1) derived by 
numerical variational means. 

The key to efficient handling of translationally invari¬ 
ant systems are transfer matrices. Their usefulness in han¬ 
dling MPS representations of infinte systems has been 
known for a while [14]. Impressive progress has been achiev¬ 
ed in describing elementary excitations within the MPS 
framework [15-18]. Haegeman et al. presented a momen¬ 
tum space ansatz [15,17] that yields very accurate results 
for the dispersion relation and is related to the calcula¬ 
tion presented in this paper. They also proved by means of 
the Lieb-Robinson bounds [19] that an excited momentum 
eigenstate of a lattice Hamiltonian can be exponentially 
well approximated by acting on the ground state with the 
momentum superposition of a local operator with finite 
support [16] if the system displays an energy gap. 

In the present paper we show how such an operator can 
be constructed using eigenvectors of an eigenvalue prob¬ 
lem arising in the dispersion calculation. We demonstrate, 
that this operator is a representation of the local creation 
operator a' by using it to compute the one-particle contri¬ 
bution to the spectral weight. To this end, we think that a 
brief presentation of the concept of matrix product states 
is required in order to explain technical details and the 
employed notation. 


2 Model and exact results 

The TFIM [21] is a common toy model for studying quan¬ 
tum magnets. In the one-dimensional ferromagnetic case 
considered here, it is given by the Hamiltonian 

u = -rY J St- JJ2 S ? S ?+ 1 - r ' J > 0 ( 2 ) 

i i 

where S z and S x are components of the standard spin- 
5 operators, T is the external field and J the coupling 
strength. We define the dimensionless parameter 


that controls the system’s behavior and in terms of which 
the Hamiltonian reads 

tt A :=-5>f-2A£SfS? +1 . (4) 

i i 

At the point A = 1 the coupling to the external field is of 
the same strength as the nearest-neighbor coupling, giving 
rise to a quantum critical point. 

The phase with A > 1 is called the Ising regime or or¬ 
dered phase since the Ising interaction is dominant. The 
ground state is twofold degenerate since its long-range or¬ 
der spontaneously breaks the Hamiltonian’s Z 2 symme¬ 
try S x —> —S x . In ID, the elementary excitations are 
(dressed) domain walls between regions of the two dif¬ 
ferent ground states. 

Conversely, for A < 1 the external field dominates the 
behavior. This phase is called the strong-field regime or 
disordered phase. Due to the strong external field, its 
ground state is unique and its excitations are (dressed) 
spin flips. The model is analytically solvable by a sequence 
of Jordan-Wigner, Fourier and Bogoliubov transforma¬ 
tions [22] as shown by Pfeuty [20] and has been studied 
extensively, see for instance Refs. [23-26]. Here we will 
recall the important known facts. 

The ground state energy per lattice site is given by the 
elliptic integral 


We test the method on the TFIM which is analytically 
solvable [20] and well understood. Inspite of its simplicity, 
it shows a variety of interesting features such as a quan¬ 
tum phase transition, ground state degeneracy and two 
different types of elementary excitations. 

The paper is structured as follows: In Section 2 the 
model and its exact solution are recalled. In Section 3 a 
short introduction to matrix product states is given, which 
may be skipped by readers familiar with this concept. Sec¬ 
tion 4 shows the results for the ground state energy and 
the dispersion followed by the derivation of the effective 
model in Section 5. In Section 6 we compute the static 
one-particle spectral weight as an application of the ef¬ 
fective model and finally conclude the paper in Section 
7. 


which displays a non-analyticality at A = 1 due to the 
singular derivative of the square root function. 

The one-particle dispersion relation uj q reads 

oj q = r yj\ + A 2 — 2Acosg (6) 

with the lattice constant set to unity. From this the exci¬ 
tation gap A is read off to be 

A = mine j q = F|1 — A| (7) 

9 

which vanishes at A = 1. 






Frederik Keim, Gotz S. Uhrig: Effective One-Dimensional Models from Matrix Product States 


3 


A quantity of great interest in connection with elemen¬ 
tary magnetic excitations is the ground state spin-spin 
correlation function defined by 

G f := (SZSf) = {fo\SSSf\iM . (8) 

where a, /? £ {x, y, z, +, —} and |^q) denotes the ground 
state. For translationally invariant gapped systems with 
local Hamiltonians in one dimension the ground state cor¬ 
relation function Eq. (8) is known to show exponential 
behavior [27] 

Gj ocexp ( 9 ) 


where £ is the correlation length. For the ID TFIM it can 
be calculated analytically [28] to be 


z = 


1 

I In A] ' 


( 10 ) 


The standard approximation £ « where v is obtained 

by fitting uj q ~ sj'A 2 + (2usin (§)) 2 to the minimum and 
maximum of the dispersion is in very good agreement with 
Eq. (10) for A > 0.2. 

Another important quantity in relating theoretical re¬ 
sults to experiment is the dynamical structure factor (DSF) 
[29] 


i _ no o 

S a ^{co,q) = d te l ^ t+q{ri ~ r ^ (S?(t)S?( 0)) 

^•>3 

(ii) 


which describes the intensity in neutron scattering exper¬ 
iments. We consider only zero temperature where the an¬ 
gular brackets denote the ground state expectation value. 
Integrating Eq. (11) over frequency yields the static struc¬ 
ture factor 

S a P(q) = |V e iq ( ri - r *US?S?) (12) 

L z ' J 


which is the Fourier transform of the ground state spin- 
spin correlation function Eq. (8). 

If for given momentum q the energy levels are well 
separated, i.e. the spectral function is a sequence of Dirac- 
(5-spikes, Eq. (11) can be written in the spectral form 

S aP {iJ,q) = Y l 8{u-E A ) s ¥(q) . (13) 

A 

The spectral weights (q) are given by projecting Eq. 
(12) onto the states with energy Ea■ Note that the energy 
Ea is understood to be defined relative to the ground state 
energy. 

Except at criticality, the TFIM has a well defined one- 
particle energy ui q . Thus evaluating Eq. (13) at Ea = 0J q 
is valid and yields the one-particle spectral weights 

S^(q) = \ |V-o)e i9W '- j) (14a) 

jd' 

= (V’o \Sq^ alli’o) {4’o\a q Sq\4’o) (14b) 


where | <j> q ) is a one-particle state with energy to q . Hamer et 
al. have given an analytic formula for the spectral weight 
in the sx-channel Sfg(q) in the disordered phase. The 
expression has been conjectured by them from high order 
series expansion [30]. In fact, its Fourier transform had 
been derived exactly by Vaidya and Tracy [31]. It reads 

s?%(q)= (i : A ;r * a < 1 ■ < 15 ) 

P w(g, A) 

Note, that there is no single-particle contribution in the 
xx-channel for A > 1 [31]. 

3 Matrix product states 

3.1 Definition 

The formalism of matrix product states (MPS) has been 
introduced in various contexts [14,32,33]. It is a way of 
denoting quantum mechanical states that is particularly 
convenient for variational calculations. It is also closely re¬ 
lated to the DMRG method [2,34,35]. This section gives 
a brief introduction to the concept. Since we are inter¬ 
ested only in translationally invariant chain models, we 
will restrict ourselves to those here. For a more detailed 
overview, we refer the reader to Ref. [2] . 

Consider a state | ip) of a system with L sites where <jj 
defines the local state at each site i 

IV>) = c <Tl> ...^ L \o 1 ,...,o L ) = 5Zc {cri }|{o-i}) ( 16 ) 

01,-..,CTL {<7;} 

where the \oi represent an ortho normal basis 
set. For simplicity, we assume that all <jj have the same 
local Hilbert space dimension d. The expansion coeffi¬ 
cients c 0 . li .., ;Crt can be interpreted as elements of a matrix 
'id w , of dimension d x d i_1 . This can be written 

as 


q,m = U m S K V Ui] (17) 

by means of the singular value decomposition (SVD). In 
Eq. (17) I/M is a d x d unitary matrix, f/W is a d x d 
real diagonal matrix that holds the singular values of 
and I/M is a d L 1 x d column orthogonal matrix, i.e., 

ytwyii] = l d . 

Now one can define the elements of the d x d L ~ l matrix 
gldyt [i] as elements of a new d 2 x d L ~ 2 matrix >f4 2 ] 


(S [ V t[11 )a ll(CT 2,..,^) “Sf 1 


ai ,o" 2 ), (cr 3 ,.. ., 0 - 3 ) 


(18) 


(with = 1,..., d), and apply the SVD again. This pro¬ 
cess can be iterated for all quantum numbers cq. In the 
end, one has 


c {<n} = (Ul 1 } ■ ■ ■ Ul L L Zl ] ■ ■ (19) 

As seen in Eq. (18), in each iteration, the quantum number 
Oi is shifted from the column index to the row index of . 
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Therefore, the matrices t/M are of dimensions d l x d l . In 
order to carry out the matrix product, one has to select the 
right block (labeled by Oj) from each f/M. In other words, 
each E/W and also the leftover can be interpreted as 
a column vector of d sub-matrices of dimension cf -1 x d l , 
which are indexed by <7,; 


U M - 



With A a ’ e C d ' lxd ' 


( 20 ) 


This quantity can also be seen as a tensor of order three 
Eventually, a single coefficient is repre¬ 

sented in the form 


<V<} = E 

ai,...,a L 

= A ai - A* 2 ■■■ A aL . 


. .. A ' 71 - -1 A ' 7L 

012 OL-l.OL £*£,1 


( 21 ) 


The index cq denotes the physical state of the correspond¬ 
ing quantum number and is therefore referred to as phys¬ 
ical index. The index 1 in A^ 1 and A aL implies that these 
objects are vectors. Finally, the entire state reads 


M = '£(A^...A^)\{a i }) . ( 22 ) 


Since each coefficient C{ ai j in Eq. (21) has the form of a 
product of L matrices, the representation Eq. (22) is called 
a matrix product state. 

By construction from the SVD, seen from the left end 
the matrices have increasing dimension: 1 x d, d x d 2 ,... 
up to i = A and are therefore distinct for different cr,. 
Since we consider ID chain models, this corresponds to 
open boundary conditions (OBC) where the position in 
the chain matters. To implement periodic boundary con¬ 
ditions, all matrices have to be of the same dimension 1 , 
because the labeling of the sites can be shifted arbitrarily. 
This is reflected in the cyclic property of the trace opera¬ 
tion, which yields the proper scalar coefficient C{ CT .} in this 
case [2] 


c {rJt} =Tr (A ai ■ A 72 ■■■A ,JL ) 

= Tr(A CT2 • • • A aL ■ A ai ) . (23) 

Since Eq. (21) is a scalar expression, applying the trace 
does not change it and Eq. (23) also holds for OBC. 

In summary, generally the (maximum) dimension of 
the A ai grows as d L ! 2 with L and may vary with the lat¬ 
tice site depending on the boundary conditions. However, 
for variational calculations, fixing all matrices to a given 
dimension D x D provides a way of truncating the Hilbert 
space which is systematic in the sense that it influences all 
bulk matrices in the same way. This D is sometimes re¬ 
ferred to as ‘bond dimension’. Then it is more convenient 

1 In general, this dimension is d L ^ 2 x d L ^ 2 . In special cases 

an exact MPS representation can be found even with 2x2 

matrices, e.g., for the ID AKLT valence bond crystal [2,14,36]. 


to have the same dimension also at the ends of the sys¬ 
tem, regardless of the boundary conditions. To this end, 
the handling of the boundary conditions is shifted to two 
auxiliary systems of dimension D located at both ends 
of the chain with states |a) and |/3). The corresponding 
matrices are vectors d“^ and bP of dimension D. Putting 
everything together results in a very general ansatz for a 
MPS 

iv>) = E E ( 24a ) 

«>/3 {o-;} 

= EE Tr(a“ t H‘ 71 • • • A aL ¥) |{tr i })|a)|0) . (24b) 

{cr;} 


Note, that in this representation the trace operation is re¬ 
dundant. However, it is still helpful in understanding the 
way matrix elements and overlaps are computed in the 
thermodynamic limit, therefore, it is kept in the notation. 
Another commonly used notation hides the boundary con¬ 
ditions in a boundary operator Q in terms of which the 
MPS ansatz reads 

li 3 ) = E Tr(QA' 71 • • • A° L ) |{cTj}) (25) 


where the trace operation is then required to make the 
coefficients scalars. 

Note that the MPS representation Eq. (22) is never 
unique. The construction starting the SVDs from the left 
side described in this section yields the so called left canon¬ 
ical form of a MPS. One could equally well start the de¬ 
composition from the right side or from both sides simul¬ 
taneously meeting somewhere in the middle of the chain. 
This would yield different matrices A ai . These canonical 
forms are very special representations because there are 
many gauge degrees of freedom generally. Between any 
two matrix sets A ai , A ai+1 one can always introduce an 
invertible matrix X such that 

c {(Ti} =Tr (A ai ■ ■ ■ A a ‘lA ai+1 ■■■A tTL ) 

= Tr (A ai ■ ■ ■ (A ai X i )(X~ 1 A ,Ti + 1 ) ■ ■ ■ A aL ) 

= Tr (A ai ■■■A' 7 'A 7 ^ 1 ■ ■ ■ A 71 - ) (26) 

which changes the adjacent matrices but leaves the coeffi¬ 
cient C { ai j unchanged. Thus, equality of two states \ifi) = 
|^ 2 ) does not imply equality of their respective MPS ma¬ 
trix sets. Therefore we understand the equality of two ma¬ 
trix sets A ai and A ai up to such a gauge transformation 
and as a shorthand meaning both sets represent the same 
state. 


3.2 Local operators 

Having defined the matrix product representation of quan¬ 
tum mechanical states, a compatible definition of oper¬ 
ators is introduced as well. Analogous to a state being 
defined by its expansion coefficients with respect to the 
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(V’olV’o) 


E EE““' t A s ‘ i* ■ • • A s ' L *b^'* a^A Sl ■ ■ ■ A SL b p {s\, 

ol', at,(3 ',/3 { Si } {sj} 

EE Tr (d aT A A * • • • A^ *b^*)Tr(a af QA Sl ■ ■ ■ A SL b?) 

<*,P Hi} 




-EE Tr \{d aT ® a Qt )(A Sl * ® A Sl ) ■ ■ ■ ( A SL * ® A SL )(b ® b B ) 

a,/3 {si} 


= Tr 


at E ASl * ® ^ I ' ’ • E A ° j * ® ^ ' E ® A Sl J b 


\sl = 1 


=:Tj 


= :T, 


=:Ti 


(29a) 

(29b) 

(29c) 

(29d) 


basis |{<Tj}), an operator O can be defined by its matrix 
elements 

({cr'}|6|{crj) = Tr (W^ ai ■ ■ ■ W a > L ) 

= W&'-WSb.a, (27) 

ai,-,aL 

where the quantities W^ 1 , iQi are tensors of order 4 and 
for given tr', cq the tT°' iCri represents a matrix. The general 
derivation and treatment of these objects is called matrix 
product operator (MPO) formalism [2]. 

Generic Hamiltonians (2) consist of terms acting only 
on a small number of lattice sites. In the TFIM one or two 
sites are involved. This simplifies the general definition 
in Eq. (27). Let O be an operator that is the identity 
everywhere except at site j, i.e., O is a single-site operator. 
Then its matrix elements with respect to two MPS are 
given by 

(0|O|^)=^Tr (ot*®ot) 

ck/3 _ V o-i J 

F a 'j* ® A aj I • • • (b* g ® b ^ 
a i a i J 

(28) 

where the matrices F ai describe the state | </>), the ma¬ 
trices A ai the state and ® denotes the Kronecker 
product. In the local Hilbert space of the single site j the 
W a j a3 are just the elements of the matrix representation 
of O, i.e., scalars. This scheme readily extends to oper¬ 
ators that are products of a finite number of single-site 
operators (see Eq. (40)). 

3.3 Thermodynamic limit (iMPS) 

Let us consider the case, where the local Hilbert spaces at 
all sites refer to locally identical spin degrees of freedom 
in a spin chain model such as the one defined in Eq. (2). 
The labels cq run everywhere over the same set of values. 


Furthermore, we assume that the Hamiltonian acting on 
these degrees of freedom is the same at each site. Then, the 
chain is translationally invariant in the thermodynamic 
limit L —>■ oo. 

Given translational invariance, it is plausible to assume 
that a uniform MPS representation exists for separable 
ground states, i.e., states that can be separated in two 
blocks by a Schmidt decomposition [14]. This means that 
all matrices A ai can be chosen the same. Such a state also 
results as a fixed point in the infinite system DMRG algo¬ 
rithm [34,35] and is called an iMPS. The uniform ground 
state matrices will be labeled A s henceforth. Next, we con¬ 
sider the norm of the ground state in Eq. (29). 

see equation (29) above. 

In Eq. (29d) we defined the boundary vectors a) := 
J2 a a“ T ® anc ] := bP* ®bP. The object T, which 
is also defined in Eq. (29d), is called transfer operator or 
transfer matrix [37] . Because the A Si are the same at each 
site the transfer matrix is also uniform: X) = T V*. 

The trace operation is redundant and used only to 
motivate the Kronecker product structure because of the 
identity Tr(H)Tr(H) = Tr(H ® B). Finally, the norm can 
be cast in the form 

<V#> =a t (T t )7T^6 , (30) 

which explains the name transfer matrix: If at some site 
b represents the right end of the chain the application of 
T transfers the this chain end by one site to the left, i.e., 
it adds the next site. 

From the definitions in Eq. (29d) it is obvious, that 
the transfer matrix T is of dimension D 2 x D 2 and the 
vectors a and b are of dimension D 2 . They can also be 
interpreted to be D x D matrices a and b by filling such 
a matrix from top to bottom and left to right with the 
vector components 

(“A 

a = I I i->- a = (ai ... a_o) (31) 

\ a D J 
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where the eq are D dimensional column vectors. In this 

n 2 

notation, the standard scalar product in C reads 

(a, b ) = Tr(a^6) . (32) 

The application of T to b or of to a is also very 
concise if a and b are denoted as matrices 

d 

T[b\=^ j A s bA s ' s (33a) 

S= 1 

d 

Tt[ a ] =^A st oT (33b) 

5 — 1 

yielding again D x D matrices. Note that, if 2d < D , this 
evaluation of these expressions is computationally more 
efficient than multiplying a D 2 x D 2 matrix to a D 2 di¬ 
mensional vector. 

Let fii be the eigenvalues, v- L the corresponding right 
eigenvectors (or eigenmatrices in the notation (31)) of T, 
and Ui the left eigenvectors, which are also the co-vectors 
of the Vi . If T is hermitian, i = Ui holds. But this is 
generally not the case. We consider the decomposition of 
b into the Vi 

b = Y,(ui,b) Vi =: ^2 

i i 

=> T[b] = ^2 PiViVi ■ (34) 

i 

By the same argument that supports the power method 
for finding eigenvalues one realizes that for very large L 
(30) implies 

= q*qu\ voPo (35) 

Mo 

where /To is the largest eigenvalue in absolute value of T 
and uq and Vq are the corresponding left and right eigen¬ 
vectors. This holds under the two conditions that the over¬ 
laps ao = (a|uo) and /?o = (6|uo) are finite and that |/to| 

is unique, i.e., there is no other eigenvalue of the same 

absolute value. 

If these two conditions are met, the explicit form of a 
and b and thus the boundary conditions they describe are 
irrelevant. From the physical point of view, this is under¬ 
stood for correlations of finite correlation length. Then, 
the behavior in the bulk of the infinite system does not 
depend on the boundary conditions. 

Note that the conclusion on the irrelevance of the bound¬ 
ary conditions does not hold for degenerate ground states 
where the boundary conditions may indeed influence the 
state in the bulk (see Appendix C). Then the exact trans¬ 
fer matrix T generically displays two or more eigenvalues 
of the same absolute value. For a physical example we re¬ 
fer to the Majumdar-Ghosh model [38-40]. In its ground 
state, two spin-4 couple into a singlet state either on the 
odd or on the even bonds which leads to a two-fold ground 
state degeneracy in the infinite system. This degeneracy is 
broken if there is a boundary: If there is a boundary, the 


realized ground state favors a singlet on the last bond in 
order to avoid a dangling spin. In the corresponding an¬ 
alytical transfer matrix we observe a two-fold degeneracy 
of \no |, i.e., the above mentioned conditions are not met. 

In the Ising phase, the TFIM also has a two-fold degen¬ 
erate ground state. As opposed to the Majumdar-Ghosh 
model, however, it does not have an exact iMPS repre¬ 
sentation at finite D. We observe that the ground state 
search produces either one ground state or the other. The 
superposition of both cannot be captured well by the MPS 
ansatz. Around each of the two ground states, the above 
stated conditions hold and |yuo I is unique so that we may 
omit the boundary vectors a and b from the notation un¬ 
less stated otherwise. In this sense, the description of the 
system reduces to computing /ro, Vq and uq and we will 
call vq and uq the boundary matrices. 

Once fio is known, A s can always be rescaled such that 
|/r 0 | = 1 and (35) stays finite for L —> oo. If /r 0 is positive, 
it can be rescaled to /To = 1, otherwise, a phase factor 
remains 2 . This scaling for /To £ M is implied in the sequel. 

Moreover, any scalar multiple of an eigenvector is also 
an eigenvector. Therefore, u o and Vq can be rescaled such 
that (aouo, /3 qVo) = 1- These rescaled eigenvectors are la¬ 
beled u and v and will be used from here on. Due to the 
gauge freedom, see Eq. (26)), one can find a gauge for A s 
such that either v or u equals the identity. Then the other 
eigenmatrix is a diagonal matrix with non-negativ eigen¬ 
values and unit trace. It corresponds to the reduced den¬ 
sity matrix appearing in DMRG. This gauge has some ad¬ 
vantages, see Appendix A for details), and is therefore the 
representation of choice. For further details of this canon¬ 
ical form of the infinite-size MPS (iMPS) see Ref. [42] 3 * . 

If not stated otherwise, we henceforth consider a non¬ 
degenerate tranfer matrix T with unique /To = 1 after 
appropriate rescaling. We consider the single-site operator 
O from Eq. (28) to be the identity with the local matrix 
representation l d . Using Eq. (30) we calculate the ground 


2 In this case, for every explicitly applied single-site operator 
(including T) the resulting matrix needs to be devided by /io to 
account for the phase factor. See Ref. [41] for a more detailed 
discussion of degenerate /r 0 . 

3 Ref. [42] uses a composite representation {F, A} of the MPS 
where F is a rank 3 tensor that lives on the sites and A is a 

diagonal matrix that lives on the bonds. It holds the Schmidt 
coefficients of the Schmidt decomposition across the bond. This 
representation can be obtained from the A B tensors by com¬ 
puting the SVD U = WAV f and setting F 3 = Ur. Here 

U is the matrix from Eq. (20) and the W s are the blocks 
of W defined analogously to the matrices A ai . The canonical 
A“ is obtained by A“ = F S A where {F, A} is the canonical 
composite representation. 
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state expectation value in the thermodynamic limit 

^(l d ) ss ^ s '* ® A s J T"-b 

s,s' J 

(36a) 

J2 A s * ® A s ^j v (36b) 

= u^Tv = (u,T[v]) = l . (36c) 

In this sense, T can also be perceived as an identity oper¬ 
ation at one site 

T[v] = A " v ^ St = 51 &a,>A a 'v A* 

s s,s' 

= 55(1 d ) ss 'A s 'vA A =: l (A ’ A) [v] = V ■ (37) 

s,s' 

The scheme in Eq. (33) extends to nontrivial operators 
straightforwardly as follows 

6 (B ’ A) [„] = ^ o ss - A s 'v R st (38a) 

s,s' 

O^BA) [ u ] = 0\ s , AA u B s ' . (38b) 


where hi is the local term of the Hamiltonian defined in 
Eq. (39). There are various methods of finding an optimal 
A s for given D. 

Eq. (41) is a highly nonlinear function in the elements 
of A s . Thus, for its minimization, one may think to resort 
to any multi-dimensional minimizer that does not rely on 
derivatives. But the convergence of them is usually slow. 
Another non-variational possibility for the ground state 
search is the imaginary time evolution in Vidal’s iTEBD 
[43] which, however, is also found to much slower than 
MPS-based iDMRG [44], 

Alternatively, one may use an iterative approach. In 
each step, only the elements of the matrcies B s at a sin¬ 
gle site are varied, all other sites are kept at a fixed A s . 
Then, the matrix B s with the lowest “local energy” eo 
is adopted everywhere as improved guess for A s and this 
process is iterated untill convergence R(| = A s is reached 
within some numerical tolerance. 

Let |i p(A s ,B s )) be the state that has A s matrices ev¬ 
erywhere except at site i = 0 where the B s matrices are 
inserted instead. This insertion breaks the uniformity of 
the state. Therefore, its energy is no longer given by a 
multiple of the expectation value of hi. Instead the full 
Hamiltonian TL = XX has to be taken into account. 
The minimization problem in terms of the elements of B s 
reads 


(V’olllV’o) = a f (T f ) * 



As an example for the application of a local MPOs let us 


consider a single term 

hi = - rS ? - JSfSf +1 (39) 

of the Hamiltonian (2). Applying the scheme in Eq. (28) 
using (38) yields 

{<t>\hi\ip) = - - J(<l>\SfSf +1 \il,) (40a) 

= -r(u,S*( F ’ A) [v]) + (40b) 

- J{u, S x (F ’ A) [S x {F ' A ) [„]]) (40c) 


where v and u are the eigenvectors of T = XX S E s * ® A s . 

This concludes the brief formal and technical review 
of the matrix product formalism. Below we turn to its 
application to the TFIM and to the construction of local 
creation and annihilation operators. 


4 Ground state energy and dispersion 

In this section, we describe one of several ways to obtain a 
uniform iMPS representation of the ground state and how 
to calculate the dispersion of a single quasi-particle in the 
system. 


4.1 Ground state search 


Starting from the ansatz (24), finding the ground state 
energy is a variational problem in the coefficients of the 
ground state matrices A s 


Eo / . (MA s )\hi\MA s )) 

— \ min- 

L "(AM (MA s Mo(A s )) 


(41) 


= {rP(A a ,B a )\^: i h i -E(A‘>))\il>{A a ,B‘>)) 

(iP{A s ,B s )\tl;(A s ,B s )} [ ’ 

where E(A S ) is the energy per site of the uniform state 
that has only A s matrices. This is also the best estimate 
for the ground state energy per site E 0 /L at each step of 
the iteration. We subtract it in order to avoid extensive 
contributions. 

Both the numerator and the denominator of Eq. (42) 
are bilinear forms in the d ■ D 2 -dimensional vector B that 
holds all the elements of the matrix set B s . Looking for 
minima in Eq. (42) means looking for roots in its deriva¬ 
tive with respect to . It is well-established that for bi¬ 
linear forms the roots of this derivative amount up to the 
generalized eigenvalue problem (EVP) 

(42) <^> B^M(H, A S )B = eB^N(A s )B (43a) 

-»■ M(H, A S )B = eN(A s )B . (43b) 

Note that the matrices M and N are both Hermitian by 
construction. For details on the ground state search algo¬ 
rithm, we refer the reader to appendix A. 

There is no rigorous proof that adopting the local min¬ 
imum Rg, that is found from the minimization at one site, 
at all sites will lower the total energy. But empirically it 
is found to be the case if the initial guess for the A s is not 
too far away from an optimal A s . Moreover, in practice we 
adopt a line-search algorithm between Rq and the former 
A s to stabilize the minimization, see appendix A.2. 

Results for the ground state energy of the TFIM are 
depicted in Fig. 1. The agreement is extremely good in 
view the low matrix dimension. There is a clear maximum 
in the deviation from the exact result, close to the location 
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Figure 1. (Color online) Upper panel: Ground state en¬ 
ergy per site Eq/L as function of A. Comparison of the ex¬ 
act result (5) to results from iMPS calculations with vari¬ 
ous bond dimensions D. The lower panel shows the deviation 
\AE\ = I So /i- -Bo,exact | ■ The critical point is located at A = 1 
and the shaded area to its right marks the Ising regime with 
two-fold degenerate ground state. 


of the phase transition. The parameter value of the largest 
deviation is found to be below the true critical values A < 
A c = 1, but quickly approaches it as D grows. 

As explained for instance in Ref. [42] , the MPS formal¬ 
ism inherently implies exponentially decaying correlations 
because of the finite, bounded amount of entanglement 
which can be represented. The length scale of the decay 
of these correlations, the correlation length, is determined 
by the second largest magnitude eigenvalue Mi of T. To 
see this, consider the application of T to a matrix b ex¬ 
panded in eigenmatrices Vi of T in Eq. (34). Assuming a 
non-degenerate spectrum of T and Mo = 1 the subleading 
term is /3i/iiTi with |mi| < 1. This term determines the 
rate at which T? [ 6 ] converges to 

T J [ 6 ] « /3 0 v 0 + n{PiVi . (44) 


Therefore, the correlation length £ 7 - captured by T is given 
as 


£r 


1 

I In Mi I 


(45) 


Figure 2 displays £t for various matrix dimensions D 
in comparison to the exact expression Eq. (10). Especially 
close to criticality a larger matrix dimension is required 
to improve the agreement. Since this is a proof-of-concept 
study, the computations were carried out on laptop com¬ 
puters and workstations. Therefore, we restricted the bond 
dimension D to low values to keep the runtime short. Bond 
dimensions of several hundred are possible, but then the 
calculations take considerably more time. 


4.2 Dispersion 

For simplicity, we focus here on the regime where the 
ground state is unique. In Appendix C the changes for 



Figure 2. (Color online) The correlation length £t as com¬ 
puted from the second largest EV of the transfer matrix, see 
Eq. (45), compared to the exact expression £ = | In A | — 1 from 
( 10 ). 


degenerate ground state are summarized. The approach is 
also described in Ref. [17]. 

If the ground state is unique, the eigenvectors B a> 0 of 
(43b) with higher local energy e a>0 describe excitations 
of the system. Let 

| ip?) = Tr (^ Sl • • • A a *~ 1 B%A a *+ 1 ■ ■ ■ A^)|{s s }) (46) 

{*.} 

be the state that has ground state matrices everywhere 
except at site j where a B s a matrix is inserted instead. By 
construction, these states are orthogonal to the ground 
state 

(*l>oWf) = («, 1 (j4 ’ Bq) M) = A*NB a = 0 . (47) 

The same holds for states with different a if the 73* are 
at the same site since 

{1>?\ ^)=BlNB p0 c5 aP (48) 

because they result from the same generalized EVP (43) 
with Hermitian matrices M and N. 

But if 23® and Bp are inserted at different sites j' ^ j 
the corresponding states will not be orthogonal. If we want 
to view the insertion A s —>• B s a as the effect of a creation 
operator we need that the excitations at different sites 
are mutually orthogonal. How can we solve this issue? We 
achieve orthogonality by resorting to the the construction 
of Wannier states known from solid state text books. One 
takes a detour via the Fourier transform since the resulting 
momentum eigenstates 

\r q ) ■■= ( 49 ) 

v 3 

are known to be orthogonal in momentum space (') (x 
5 qq ' because they refer to different eigenvalues under dis¬ 
crete translations. The restriction of all momenta to the 
first Brillouin zone is implied. 
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Exploiting translational invariance, the overlap of two the dispersion relation can be found by a second varia- 
momentum eigenstates can be computed as tional calculation in the orthogonal complement of the 

i ground state, i.e., in the sub space spanned by the |i/;“ >0 ) 



(50a) 


3,3' 


z . (<l> q W-E 0 ) |0,) 

u a < mm —-—;—;— - — 


(50b) 

{4 > q\ ( f > q) 

jJ' 


d-D 2 — l 


(50c) 

:= E ■ 

a.=l 

i a P 

q 

(50d) 

This leads to another generalized EVP 


(53a) 

(53b) 


where Eq. (50) defines the matrix N q which is the metric 
tensor of the states | ip q ). As seen in Eqs. (50), the nor¬ 
malization factor ) always cancels out in the computation 
of an expectation value, norm or overlap. Thus the limit 
L —> oo does not pose any numerical problems. As seen in 
the sequence of equalities (50), translational invariance al¬ 
lows us to assume that the ket-side matrices B s are always 
placed at site 0 . 

The infinite sums over real space indices in the above 
equations may be seen as insurmountable problem. But 
this is not the case because the infinite sums converge 
exponentially. To elucidate this point we look at the fol¬ 
lowing limits. Any D x D matrix to (assuming /io = 1) 
can be decomposed according to Eq. (34). Applying the 
transfer matrix T j-times yields 


T j [m\ = E (ui,m)T jv i 

i. 


II 

-M- 

e 

tT 

lS». 

(51a) 

with \ni\ <1 for i > 0 


=> lim T° [to] = (it, to) • v 

j-to o 

(51b) 

lim = ( to , v) ■ u 

i->oo 

(51c) 

where the convergence to these limits is exponential in j 
governed by the second largest absolute value \fii\ of the 
eigenvalues of T. 

The overlap (V , “|V’o) § oes to zer0 ^ or l &r g e j 


(52a) 

w {l (Ba ' A) [u],(u,l {A ’ Bp) [«])•«) 

(52b) 

= (V’o|V , o)(EI^’o) 

(52c) 

= 0 

(52d) 


where the approximation in (52b) refers to the exponential 
convergence established in Eqs. (51). In the last line we 
used the local orthogonality (V’oll/'o ) = 0. The exponential 
convergence to zero justifies to trunate the Fourier series 
after a finite number ) max of terms. 

If the ground state search is converged, B a=0 is the 
set of ground state matrices, i.e., |^“ =0 ) = |^o) V i. If the 
ground state matrices are established with sufficient nu¬ 
merical accuracy, all |V'“ >0 ) are orthogonal to |i/>o). Thus, 


H q V q = UJqNqVq (54) 

where N q is the matrix defined in Eq. (50d) and H q is 
defined analogously 


Hf := (1>Z\(H-Eo)\1>t) 

= Y j e iqj hf i 




(55a) 

(55b) 

(55c) 


The lowest eigenvalue is the best estimate for the one- 
particle dispersion at given q. 

Haegeman et al. [15] observed that there is always a 
number of choices £?* such that \i/jq) = 0 V q due to the 
gauge degrees of freedom stated in Eq. (26) combined with 
translational invariance. Because of the associativity of the 
matrix product for any A' £ C nx " we have 

IVf-i) = E Tr (- • • ■ • -)|{*}> = 

{«<} 

\^) = J2M---A s ^(XA^)...)\{ Sl }) . (56) 

RI 


Note that we allow here for the more general case where 
the ground state matrices are different to the left and to 
the right of the inserted gauge matrix X. This includes the 
possibility of excitations of domain wall character where 
one switches between degenerate ground states. 

Let us define B s := e lq A s X — XA S implying 



(57a) 

V 3 

(57b) 

V 3 

(57c) 

0 

(57d) 


because a phase factor of e iq translates to a shift of the 
states in real space by one site to the left under the Fourier 
transformation. The matrix X has D 2 parameters. Thus, 
for q ^ 0 or for q = 0 and A s ^ A s , the dimension of the 
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space spanned by \ij} q ) is reduced by D 2 . In other words, 
there are D 2 “zero modes”. 

For q = 0 and A s = A s , the choice of X = 1 results 
in B s = 0 which makes \ipj) the null vector. Thus the 
number of linearly independent zero modes is reduced to 
D 2 — 1 in this case. Therefore, the metric tensor N q has 
a D 2 or ( D 2 — 1) dimensional null space. In order to take 
this into account H q needs to be projected onto the non¬ 
zero eigenspace for computing uj q . Within this non-zero 
eigenspace, the dispersion is found by solving the standard 
EVP 


y/D^ \*H q V'yfD^ 1 v , q =u q v , q (58) 


where the diagonal matrix Pj\r holds the non-zero eigen¬ 
values of N q and V' the corresponding eigenvectors. The 
original v q from Eq. (54) can be obtained as v q = V'v'. 

The computation of the matrix elements in Eq. (55a) 
is the most time consuming part of the calculation because 
the complete Hamiltonian acts on all lattice sites and the 
Fourier coefficients have to be computed for many values 
of j. But by the same argument as in Eq. (52), the contri¬ 
butions converge exponentially to zero if \j\ 1. Let for 

instance j <C 0, i > 0. Then 

h$ = (59a) 

(59b) 

« (l (B “' A) [(l ^-^[u]^) ■u],T < - 1 [fc i [i;]]) (59c) 

= (V’olV’o)(l (B “- A) [u],T® _1 [fti[t;]]) = 0 . (59d) 


where the vanishing of the last expression holds in the 
limit j oo. 

For |i| 1 we obtain similarly 

hfi = (i( B “>- A )[Ttl J ' _ 1 l[ 1 (- A,B P) [ u]]],T i_ 1 [ /i i[ u]]) 

- 

- 

= (V’jlV’o)^ - ^(V’jlV’o) = 0 for j -too (60) 

where j < 0 is assumed for simplicity. 

Figures 3 through 5 show the dispersion u} q for various 
parameter values. At A = 0.8 and A = 1.2, see Figs. 3 and 
4, the agreement is very good, both in the strong field 
and in the Ising regimes, because the system is placed 
not too close to criticality. The nice agreement illustrates 
that ground state degeneracy is handled very well by the 
advocated method. 

Directly at the quantum critical point, see Fig. 5, the 
closing of the gap is difficult to capture numerically. The 




Figure 3. (Color online) Dispersion for A = 0.8 in the strong- 
field regime. Comparison of the exact result (6) to results from 
iMPS calculations with various matrix dimensions D. 



q / it 

Figure 4. (Color online) Dispersion for A = 1.2 in the Ising 
regime. Comparison of the exact result (6) to results from iMPS 
calculations with various matrix dimensions D. 


reason is the diverging correlation length £, see (10). The 
amount of entanglement that can be described by an MPS 
is bounded by the matrix dimension D. Thus, no finite- 
dimensional MPS can completely describe a state with 
diverging correlation length. 

The inset in Fig. 5 depicts the gap A as function of A 
in the vicinity of A = 1. Note that the gap values are as 
low as 10” 6 P to 10 _5 P in spite of the limited bond dimen¬ 
sion. The occurrence of a rather sharp minimum indicates 
a possible phase transition. Note that this criterion is in¬ 
dependent of a comparison to the exact result and allows 
one to estimate the corresponding critical parameter value 
as well. 


5 Effective model 

As mentioned in Sect. 2, in the strong-field limit A —> 0 the 
elementary excitations are flips of single spins from the po¬ 
larized ground state. If the Ising interaction is switched on 
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q Ik 

Figure 5. (Color online) Dispersion at A = 1, i.e., at the 
quantum critical point. Comparison of the exact result (6) to 
results from iMPS calculations with various matrix dimensions 
D. The inset focuses on the gap as function of A around the 
critical point. 


by a finite A, such a flipped spin acquires a virtual dress¬ 
ing which corresponds to a polarization cloud around it. 
This means that the elementary excitations are no longer 
strictly local, but “smeared out” over a certain region on 
the chain. This concept has been the basis of the CUTs in 
real space representation [5,6,8,9]. The spatial extension 
of the polarization cloud, i.e., of the smeared out region, is 
governed by the correlation length. This can also be seen 
from Eq. (4) in Ref. [16]. 

In order to make progress in deriving effective mod¬ 
els (1) in terms of the elementary excitations we want 
to establish the key ingredients of second quantization, 
namely the creation and annihilation operator of an ele¬ 
mentary excitation. Thus, it is our objective in this sec¬ 
tion to explicitly derive a local creation operator acting on 
the ground state. If we know the ground state (or a very 
good numerical representation thereof) and we are able 
to characterize the local excited states we can follow the 
route advocated previously [6] to determine the effective 
model on the bilinear level. 

More work will be required for the determination of 
decay terms (lb) and two-particle interactions (lc). To 
determine them, states with two quasi-particles at sites i 
and j must be properly defined. This requires that they 
are normalized and two such states (i,j) and ( i',j') are 
orthogonal if i ^ i' or j ^ j'. Moreover, such two-particle 
states must fit to the one-particle states in the sense that 
they decompose into the one-particle states for |* — j| —> 
oo. These issues set the roadmap of research, but they are 
beyond the scope of the present article. 

In order to construct a local creation operator, we con¬ 
sider the eigenvector v q of Eq. (54) that belongs to the 
lowest eigenvalue which defines the dispersion u q . Its com¬ 
ponents v q describe how the states |i/>“) are linearly com¬ 
bined to form an elementary excited state that satisfies 


{<t>q\(^l - E 0 )\(j> q ) = u q (61a) 

\<l>q) = V q\^q) = ^IV’o) ■ (61b) 

a 

This means that | cf> q ) can be interpreted as a state, in 
which one quasi-particle of momentum q has been cre¬ 
ated. Taking the inverse Fourier transform of Eq. (61b) 
one obtains an expression for the action of a local cre¬ 
ation operator a\ on the ground state 

a l I'V>o) = v i I ^ i+j ) ( 62a ) 

j,oi 

with vf := ^-J 2 v 9 eiqi ■ ( 62 b) 

9 

This equation is the key element in advancing towards 
effective models via MPS representations. 

In the thermodynamic limit q is a continuous variable 
and the sum in Eq. (62b) becomes the integral over the 
Brillouin zone 

< = < 63 > 

Although numerical integration always comes down to sum¬ 
mation at some point, the continuous representation is ad¬ 
vantageous for adaptive algorithms because Eq. (54) can 
be evaluated at arbitrary values of q. See Appendix B for 
comments and technical details on handling v q . 

Taking the sum over a first in Eq. (62a) simplifies the 
numerical computation. Hence, we define a single matrix 
set 

q:='E v i B a ( 64 ) 

ex 

to be inserted at distance j from the center site i of the 
particle created by a\. In this description the particle is 
represented by a number of matrices {C!j} as follows 

«!l^o)= I ${Cj)i+j) (65) 

3 ~~ 3 max 

where \i/i(Cj)i + j) is a state analogous to |^>“) that has A s 
matrices everywhere and C? inserted at site (i + j). 

To quantify the degree of localization of the excita¬ 
tions, we study the squared norm of the vectors v 7 

Vj:=IM 2 = EW- ( 66 ) 

ex 

Figure 6 shows V) for various values of the bond dimension 
D and compares their dependence on j to the decay of the 
correlation function Gj. Clearly, the distributed (smeared 
out) contributions to the quasi-particle decay exponen¬ 
tially with the distance j from the center site. This agrees 
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Figure 6. (Color online) The quantity Vj defined in Eq. (66) 
for A = 0.9 and various matrix dimensions D. The solid red 
line shows the function 0.2 • exp(—j/£) where £ is the exact 
correlation length from (10). Vj displays exponential decay on 
a length scale that is always smaller than the correlation length 
£■ 

with the findings in Ref. [16]. With increasing matrix di¬ 
mension the decay becomes slower and approaches the de¬ 
cay of the correlation function. This is is consistent with 
the finding in Fig. 2 illustrating that larger D allows one 
to capture longer correlations. For the numerics, it is very 
advantageous that the decay of Vj is always even faster 
than the decay of the correlation function defined by the 
correlation length £, because this fact implies that the rep¬ 
resentation can be truncated after a fairly small number 
of sites |j| < j max . 


6 Spectral weight 

To illustrate the validity of the creation operator defined in 
(61b) we compute the spectral weight Sf .For a = j3 = x 
Eq. (14b) becomes 

STp (q) = (i’o\Sq^a] ] \ijj 0 )(ijjo\a q Sq\ipo) 



where m q is defined by m q = {ipo\a q S q \4>o). Inserting the 
definition of a] (62) and the Fourier transform of S q we 
obtain 

m q = \ Y, <V^e-^“|^ 0 > (68a) 

= ^ v T eiqri ^i\SoHo) (68b) 

i,a 

where the matrix elements I^q |?/>o) can be computed 
in analogy to the single-site operator in Eq. (40). 

Figures 7 and 8 depict the spectral weight in compari¬ 
son to the analytical result Eq. (15) for various values of A 
and D. For smaller values of A, see Fig. 7, well away from 



0 0.05 0.1 0.15 0.2 0.25 0.3 


q/it 

Figure 7. (Color online) Upper panel: The spectral weight 
Sfp (q) for A = 0.5 and various matrix dimensions D. Lower 
panel: The deviation of the iMPS results from Hamer’s formula 
Eq. (15). The plot interval [0,0.3] is chosen to emphasize the 
deviation for small values of q where Sfp ( q ) has its maximum. 



q /it 


Figure 8. (Color online) Upper panel: The spectral weight 
Sfp (q) for A = 0.99 and various matrix dimensions D. Lower 
panel: The deviation of the iMPS results from Hamer’s formula 
Eq. (15). The plot interval [0,0.3] is chosen to emphasize the 
deviation for small values of q where Sip ( q ) has its maximum. 

the critical point A = 1, the agreement is very good for all 
values of q. Still, larger values of D imply an even better 
agreement. For a value of A closer to the critical point, the 
agreement is still good, see Fig. 8, in view of the small val¬ 
ues of D. But in particular close to the almost diverging 
correlation at q = 0, larger values of D are indispensable 
to capture the correct correlations. 


7 Conclusions 

The objective of the present paper has been to sketch the 
roadmap to a derivation of effective one-dimensional mod¬ 
els by a numerical variational approach. In particular, we 
have explicitly shown how the first step works, i.e., the sys- 
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tematic construction of a local creation operator acting on 
the ground state. Thereby, bilinear terms in the Hamilton 
operator such as the dispersion can be determined, cf. the 
generic Hamiltonian Eq. (1). 

We have shown how the matrix product state (MPS) 
formalism can be used to derive effective models in terms 
of quasi-particles from microscopic local spin model Hamil¬ 
tonians. Based on transfer matrices, MPS work efficiently 
in the thermodynamic limit (iMPS). Starting point of the 
MPS is the accurate determination of a MPS representa¬ 
tion of the ground state. This defines the vacuum of exci¬ 
tations similar to the reference state in continuous unitary 
transformations [6]. 

A side product of our ground state search algorithm 
are eigenmatrices with higher local energies. We have shown 
how this side product can be exploited to construct the 
elementary excited states. Such constructions work very 
well for unique and for degenerate ground states. In the 
latter case the elementary excitations generically are do¬ 
main walls between the degenerate ground states. 

We derived an expression for the action of a local quasi¬ 
particle creation operator on the ground state. These quasi¬ 
particles are no longer completely local, but they are found 
to be “smeared out” but localized around one lattice site 
similar to Wannier states for the band electrons. The ap¬ 
proach is illustrated and tested for the excitations in the 
transverse-field Ising model in one dimension in the dis¬ 
ordered strong-field phase as well as in the ordered Ising 
phase. In the strong-field phase the elementary excitations 
are spin flips while they are domain walls in the Ising 
phase. 

It turns out that the quasi-particles are exponentially 
localized on a length scale always smaller than the correla¬ 
tion length £. In this way, the numerical representation of 
the elementary excitations is well controlled. Using this 
definition, the one-particle contribution to the spectral 
weight in the xx-channel has been computed. The very 
good agreement with Hamer’s formula [30] confirms his 
conjecture and strongly corroborates the validity of our 
approach. 

What are the next steps on the roadmap to effective 
models from variational approaches? In order to be able to 
determine the parts of the Hamiltonian (1) which describe 
the decay of quasi-particles (lb) or the interaction of a pair 
of them (lc) we need to extend the definition of single 
particle states to states with two-particles. The key issues 
are a proper orthogonalization of states with excitations 
at different sites. Furthermore, it must be ensured that 
the two-particle state of two very distant quasi-particles 
equals the state obtained from the successive application 
of the creation operator defined from single-particle states. 
These issues are beyond the scope of the present article, 
but represent future research. The ultimate aim is to be 
able to write down effective models in second quantization 
in terms of the elementary excitations. 

An interesting step towards this aim has been accom¬ 
plished very recently by the variational construction of 
scattering states of two elementary excitations [18]. But 


so far the explicit construction of the effective model has 
not been realized. 

A longer-term vision consists in the generalization of 
the presented approach to higher dimensions by passing 
from matrix product states to projected entangled pair 
states. The conceptual issues and their solutions, for in¬ 
stance the construction of Wannier type of local excita¬ 
tions, are the same in higher dimensions. But the numer¬ 
ical handling is less efficient than in one dimension where 
the thermodynamic limit is easily built-in by transfer ma¬ 
trices. 

In summary, we are convinced that the construction 
of effective models via numerical variational approaches 
constitutes an interesting and promising route to capture 
the physics of strongly correlated systems. 

We gratefully acknowledge the financial support of the Helmholtz 
virtual institute “New States of Matter and Their Excitations”. 
We also thank B. Fauseweh and N.A. Drescher for many help¬ 
ful discussions. 


A Ground state search algorithm 

A.l Local fixed point iteration 

This appendix contains a more detailed discussion of the 
ground state search algorithm. The starting point for a 
given guess for A s is computing fj, o, v and u. For reasons of 
both computational efficiency and algorithmic simplicity, 
rescaling is done such that /r 0 = 1 and {u,v) = 1. 

Next, we consider the normalization constraint in Eq. 
(42), i.e., the division by the norm of the state \ip(A a , B s )) 
that has A s matrices on all sites but one where the B s are 
inserted instead 

^(A a ,B a )\^(A a ,B a )) = (n,l (B ’ B) M) ± (u,v) . (A.l) 

The reference site where A s is replaced by B s is labeled 
i = 0 which marks the center of the chain. The norm 
in Eq. (A.l) is a bilinear form in the coefficients of B s , 
interpreted as a single d ■ H 2 -dimensional vector B. To see 
this, we recall how it is computed 

( W ,l( s ' B )[u]) = Tr(«t l^ B ’ B )[x]) (A.2a) 

^5 ss .B s 'vB s ^ I (A.2b) 

s,s' J 

= u a'pfiss'Bpp,vpi a (B ^) aa i 

a' s,s' (3,(3',a 

(A.2c) 

= yi yi Bj a ,(s SS 'ui ) ,,p v p' a )Bpp l . 

s,s' a,a',(3,(3' 

(A.2d) 

If one takes the three-tuples (s, a , a') =: r and (s', /3, /?') = 
t' as single indices running from 1 to d ■ D 2 , the matri¬ 
ces B s correspond to a vector B; we call this step ‘vec- 
torization’. In this notation, the norm (A.2) simply is a 
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vector-matrix-vector product 


(A.6) 


N tt ,B t , = B^NB (A.3a) 

T,T / 

with N tt i = N*? at p,p = Sss'U^pVp'a (A.3b) 

and B t = B s a p . (A.3c) 

If the vectorization B of B s is performed by concatenating 
all columns of B s into a single column (cf. Eq. (31)), a 
short calculation shows, that N is given by the d-D 2 x d-D 2 
matrix 


N := (g> v T <g) v) 


(A.4) 


where v T is the transpose of v. 


^(A s ,B s ) 0 \(n-E 0 )mA s ,B s ) 0 ) 

= £ (m s ,B °) 0 | (fh-f) mA\B°)o) 

i ——oo ' ' 

= B f (mM - N* 

%——OO ' 

+ B f + M [S -' S Z ] - B 

+ (m [ s ° i + m [s o s ri j B 

+ f>t /V^-IV^B 

2=1 ' ' 

= B*M{H- E 0 )B 



(A.7a) 
(A.7b) 


(A.7c) 
(A.7d) 
(A.7e) 


In a next step we study the application of the single¬ 
site operator S'* as it occurs for instance in Eq. (40b) . The 
algebraic structure is the same as in Eq. (A.2), but the 
identity operation is replaced by S z . Thus, its expectation 
value with respect to \ip{A s , B s )) is also a bilinear form in 
B given by 

( U ,5^ s ' s )[u]) = ^ £ B s a l(S; s ,ul 0 v p/a )B^, 

s,s' a,Oi',f3,(3' 

(A.5a) 

= B^M^B (A.5b) 

M PS1 := S z <g> v T <g> (A.5c) 

where S' 2 is the local representation of the spin-^ 2 -operator 
S' 2 = \(r z i cr 2 being the 2 -Pauli matrix. 

The Ising interaction in Eq. (40c) is a bit trickier be¬ 
cause two sites are involved. Let them be the sites i = 0 
and i + 1 = 1. Then the bilinear form reads 


where the matrices M^ hi 1 are given by 

Id 8 v T 8 
Id 8 v T 8 

with 5 = T’ l ~ x [hi[v\[ 
and u = T^~ 2 [h\[u]] 

and N is the matrix defined in Eq. (A.4). 

Note, that in Eq. (A. 7c) one summand — N^ is omit¬ 
ted. The lowest eigenvalue po then converges to the vari¬ 
ational ground state energy per lattice site rather then to 
zero. As a rule of thumb in numerics it is better to search 
for finite values than zero, especially when dealing with 
relative errors. 

As shown in Eq. (57), the boundary contributions quick¬ 
ly converge to zero for i 0 so that the sums can be 
truncated after a finite number of sites, i.e., comprising a 
finite tractable number of terms. 

Eq. (42) can be recast into the form 


if * > 0 
if * < 0 


(A.8a) 


(A.8b) 
(A.8c) 


(u,S x( - B ' B '>[S x( - A ’ A '>[v]]) = B^M^ S * ] B (A.6a) 

M [s ° s * 1 := S x 8 v T 8 ■u t (A.6b) 

v :=S x{A ' A) [v\ . (A.6c) 

One realizes in Eqs. (A.5c) and (A.6b) that the structure 
of these expectation values is always the same. It is the 
Kronecker product of the local operator matrix at site 
i = 0 where the B s matrices are inserted and the Kro¬ 
necker product of the boundary matrices u ( u ) and v (v). 
If there is only the operator at site 0 (Eq. (A.5c)), the 
boundary matrices u and v are used without modifica¬ 
tions. If there are more operators involved, see Eq. (A.6), 
the boundary matrices are modified by applying these op¬ 
erators to them. We denote this by the tilde symbol, see 
Eq. (A.6c). 

Finally, the whole numerator of Eq. (42) is given by 
the term from Eq. (A.5) and a sum of terms of the form 


B^M(H, A S )B - ^-B^N(A S )B = 0 . (A.9) 

_Zv 

The minimization of Eq amounts up to finding roots of the 
derivative with respect to Bb This yields the generalized 
EVP in Eq. (43b). 

The complete algorithm runs as follows: 

1. Start with an initial guess for A s . 

2. Generate the matrices M and N. 

3. Solve the generalized EVP Eq. (43b). 

4. Break if Bg = A s within a given tolerance. 

A local minimum of the ground state energy is a fixed 
point of the iteration which satisfies this condition. 

5. The eigenvector B 0 with lowest energy eo is chosen as 
new guess for A s . Go to step 2. 

If no better initial guess is available, start with a random 
matrix set in step 1. In cases where multiple values of a 
system parameter are to be investigated, the converged 
solution for a nearby value generically constitues a good 
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initial guess. Since there is no way to determine if a solu¬ 
tion is also a global minimum, the algorithm is terminated 
if a fixed point Bq = A s is reached within numerical tol¬ 
erance. 

In contrast, conventional MPS-based iDMRG, see for 
instance Ref. [44]) follows the spirit of White’s original 
DMRG for infinite systems [1]. In each iteration, the sys¬ 
tem is incremented by adding one unit cell in the center 
of the chain, thereby increasing the bond dimension lo¬ 
cally which is then truncated to its original value using 
the most relevant part of the density matrix. Convergence 
is reached when the matrices A* +1 obtained for the added 
unit cell are the same as the A ® obtained in the previous 
iteration within some preset tolerance. 


As can be seen from Eq. (40), the ground state energy 
is a highly nonlinear function of the coefficients of A s . Any 
multi-dimensional minimizer, that does not rely on deriva¬ 
tives, can be used to find a minimum, e.g., the method of 
conjugate directions or simulated annealing. On the one 
hand, our experience shows that these routines converge 
at a much slower rate than the algorithm described above. 
On the other hand, however, the fixed point iteration may 
fail to converge if the initial guess is too far away from an 
optimal solution. Therefore, we actually use a hybrid al¬ 
gorithm. First, a couple of iterations of simulated anneal¬ 
ing are preformed, yielding a good inital guess. Then this 
guess is used for the above given algorithm which breaks 
if Bq* = A s within numerical tolerance. 


A.2 Fine tuning 

Figure 2 shows that close to criticality the results improve 
visibly with growing matrix dimension D. This is so be¬ 
cause a larger bond dimension allows more correlations 
to be represented. At or close to the strong-field limit, 
however, a large bond dimension may actually be disad¬ 
vantageous due to a certain lack of entanglement. The 
matrices A s are too large to encode the small amount of 
entanglement in the system. This makes itself felt in the 
norm matrix N in Eq. (43b) becoming singular or having 
very small eigenvalues in magnitude. Solving the general¬ 
ized EVP involves division by the eigenvalues of N so that 
the generalized EVP is ill-defined if N is singular or close 
to it. 

But Eq. (A.4) shows that in the gauge where v = 1^, 
N is diagonal and holds d ■ D copies of u. This matrix u 
is the reduced density matrix of the left subsystem if the 
system is split into a left and a right part which is traced 
out, as is done in DMRG. Therefore, omitting the vectors 
corresponding to small eigenvalues of u is a systematically 
controlled way to focus on the relevant subspace. The null 
space of N is projected out which also avoids numerical 
instabilities, cf. Eq. (58). 

Let D' be the dimension kept in the truncated density 
matrix u. Projecting out the null space of N results in 
d ■ D ■ D' eigenvectors of Eq. (43b) instead of d ■ D 2 . This 
is also efficient in the subsequent calculations by speeding 
up the dispersion calculation since the initial dimension of 
H q and N q is reduced to d ■ D ■ D'. 

A last aspect in the ground state optimization concerns 
the iteration Bq —> A s . As mentioned in the main text, it 
is not at all clear whether taking the matrix set Bq found 
for a single site at all sites indeed improves the ground 
state. But we can ensure that the variational ground state 
energy is reduced in each iterative step by performing a 
linear search using the ansatz 

Eq = min E(cos(x)A s + sui(x)Bq) (A.10a) 

x£(0,7r/2) 

-t B if = cos(a; m i n )A s + sin(a: min )Ro (A.IOb) 

which interpolates between A s and Bq . The one-dimensional 
minimization of E{x) as function of x is numerically ro¬ 
bust. 


B Creation operator 


This short appendix contains general technical details of 
the computation of the representation of the creation op¬ 
erators aj and aj. 

The generalized EVP Eq. (54) is solved for each mo¬ 
mentum value q independently. Therefore, an arbitrary 
phase may always occur between the eigenvectors v q and 
v q +Aq where Aq is the sampling interval in g-space. In or¬ 
der for the Fourier transform Eq. (62b) back to real space 
to yield well-localized components each component v q 
must be a smooth, 27r-periodic function in q. 

In order to ensure this smoothness we employ a two 
step process. First, we fix the phase between adjacent vec¬ 
tors only separated by Aq to zero by setting 




with <P q := 


V q-Aq V q 

\Vq-Aq\W\Vq 


(B.l) 


But this choice of vanishing phase is still somewhat ar¬ 
bitrary. More generally, a phase of the order of Aq could 
occur between adjacent vectors. 

In practice, we check for the phase between v q= _ T and 
v q -„ after the above smoothing process. It should vanish 
due to 27r-periodicity, but this may not be the case. To 
restore 2-7r-periodicity in a second step, the accumulated 
phase between v q -_ n and v qz=7r is distributed evenly over 
the whole Brillouin zone. This procedure results in the 
very fast decaying quasi-particle representation presented 
in Fig. 6. 


C Ground state degeneracy 

As mentioned above, ground state degeneracy is reflected 
in the spectrum of the transfer matrix T. If an exact iMPS 
prepresentation of the ground state exists at finite D , this 
results in a degeneracy of the largest absolute value of 
the eigenvalue /io- This is for instance the case for the 
Majumdar-Ghosh model [38-40] that has an exact ground 
state iMPS representation at D = 3. If no exact iMPS 
exists for finite D , as is the case for the TFIM, one still 
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observes that 

1 1 for D —> oo , (C.l) 

lMi| 

i.e., there is an asymptotic degeneracy. 

However, in some cases analytical considerations may 
help. For the TFIM in the Ising phase we know that the 
ground state |'0 q") with magnetization in positive S x direc¬ 
tion can be transformed into the ground state \i/>q) with 
magnetization in negative S x direction by a 7r-rotation 
about the 2 -axis 

^-)=e™ s *\^) ■ (C.2) 

This rotation is a non-local operation, but it is the same 
for all sites 


C.l Dispersion calculation with degenerate ground 
state 

As mentioned in Sect. 2, the elementary excitations in the 
strong-field and in the Ising phase are qualitativly dif¬ 
ferent. In the strong-field regime, they consist of a local 
perturbation of the otherwise uniform ground state as de¬ 
scribed by the ansatz in Eq. (46). In the Ising phase, the el¬ 
ementary excitations are domain walls separating regions 
of different ground state, i.e., the excitations are non-local. 
The domain wall character requires a modified ansatz 

IO = X! Tr (^ Sl • • • A 8 *- 1 B 8 ^ A 8i+1 ■ ■ • A Si )|{sJ) 

(C.6) 


e™ sz \4’+) = Tr 


Recall that e l '* s * = cos(7r/2)l + *sin(7r/2)cr 2 = io z 
since S z = \a z where a z is the 2 -Pauli matrix and er z2 = 
1. The phase factor i is a special case of the gauge trans¬ 
formation Eq. (26) and can be dropped. In the iMPS rep¬ 
resentation the spin rotation results in a relative sign be¬ 
tween the two matrices 

I-> h/> 0 -> :{A\-A 2 }. (C.4) 

Note, that this representation of \4>o) is not canonical any¬ 
more. But it can be made canonical by the algorithm men¬ 
tioned in the main text and presented in Ref. [42]. The 
resulting density matrix u turns out to the be same as for 
IV’o")- 

The two ground states can be distinguished by the sign 
of the magnetization in /S^-direction 

M X = (^\S X \^) (C.5) 

which serves as the order parameter in the Ising regime. 
The ground state search algorithm produces either one or 
the other realization, not a superposition of both. This 
is due to the fact, that an iMPS representation strongly 
favors pure and finitely correlated states. Obviously, this 
is the case for either state 1^) but not for their super¬ 
position. Near the degeneracy of /io one eigenvalue domi¬ 
nates numerically and the algorithm converges to the cor¬ 
responding eigenvector as fixed point. Which state will 
finally be selected depends on the initial guess for A s . 

For general models, however, it may not be clear, how 
degenerate ground states are connected, i.e., if there is an 
analytically applicable transformation such as Eq. (C.2). 
But the occurrence of degenerate transfer matrices is a 
strong indicator for degenerate ground states. The com¬ 
parison of the expectation values of possible order param¬ 
eters for different ground state solutions may help to dis¬ 
tinguish them. 


A^j x... 

£ {e^ sz ) SLS , L A 8 '- 

\ s L’ s l 


!{*}> • 




where A s describes an alternative ground state. To obtain 
the appropriate eigenmatrices B a a different kind of gen¬ 
eralized EVP has to be solved once the two ground states 
are known 


M(A S ,A S ,H)B = eN{A s , A S )B (C.7) 

where the matrices M and N are built in the same way 
as M and N from Eq. (43b), see Appendix A for de¬ 
tails. But there is one important difference. Instead of the 
eigenmatrix v of T = Y s A s * <g> A s the eigenmatrix v of 
T = Y s A s * (g> A s has to be used. In this case, the eigen¬ 
matrix B a= o no longer is one of the ground state matrix 
sets A s or A s . 

For domain wall excitations it is not directly evident, 
that the \ipf) and their Fourier transforms \ip q ) are or¬ 
thogonal to the ground states although this still holds. 
This is simply due to the fact that in the thermodynamic 
limit different ground states are orthogonal. Therefore, 
any excited state that contains a domain wall and thus 
regions of both ground states is orthogonal to l^) for 
the infinite system. 

To see this in the iMPS representation we consider the 
overlap of the two ground states 

=u*Ji qV (C.8) 

where flo is the largest eigenvalue in absolute value and 
u, v are the corresponding left and right eigenvectors of 
T = A s * ® A s . Since \i/j + ) ^ |V ,_ )> Mo 7^ Mo holds. 
It turns out, that \p,o\ < 1. Therefore, the overlap (C.8) 
tends to zero as L —> oo. This implies that the ground 
states are orthogonal in the thermodynamic limit. 

Thus, the dispersion calculation is conceptually the 
same for degenerate ground states as for the non-degenerate 
case. Only the dimension of the matrices H q and N q in Eq. 
(54) may be increased by one because B a — o may also rep¬ 
resent an excitation. In the computation of the matrices 
H q and N q one has to account for the different ground 
states. For instance, the overlap of two states as defined 
in Eq. (C.6) is given by (j < 0) 

m\<) = (c . 9) 
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